      Program Fast_Fourier_Transform

C     Programm zur Berechnung der Potenzen der Einheitswurzeln

C     N  Problemgroesse; 2er Potenz

C     N2 Anzahl der zu berechnenden Potenzen; N2 = N/2
C     Direkt Fouriersynthese bzw. Fourieranalyse (1/-1)
C     W(N) Feld zur Aufnahme der Potenzen der Einheitswurzeln


      integer N,N2,k

      complex W(:)
cmf$  layout W(:serial)
      complex C(:)
      real F(:), maxf

      real Direkt

      print * ,'Input k : (N = 2**k)'
      read * ,k
      N = 2**k
      N2 = 2**(k-1)
     
      allocate(W(0:N2-1))
      allocate(C(0:N-1))
      allocate(F(0:N-1))
c
c     random initialization
c
      call cmf_random (f)
c
!HPF$ INDEPENDENT, LOCAL_ACCESS
      do i=0,N-1
        c(i) = f(i)
      end do
c
      Direkt = -1
      call Presort (N,k,C)
      call Roots (W,N2,N,Direkt)
      call Transform(C,W,N,N2,k)
c
c     inverse fft
      Direkt = 1
      call Presort (N,k,C)
      call Roots (W,N2,N,Direkt)
      call Transform(C,W,N,N2,k)
!HPF$ INDEPENDENT, LOCAL_ACCESS
      do i=0,n-1 
        f(i) = abs (c(i)/n - f(i))
      end do
      maxf = 0.0
!HPF$ INDEPENDENT, LOCAL_ACCESS
      do i = 0, N-1
         reduce (maxval, maxf, f(i))
      end do
      print *,'Max = ', maxf
      deallocate(F)
      deallocate(C)
      deallocate(W)

      end

C     *******************************************************

C     -------------------------------------------------------
C     ------------ Umordnung der Speicherplaetze ------------
C     ------------     nach Cooley und Tukey     ------------
C     -------------------------------------------------------

      subroutine Presort (N,k,C)
      complex C(0:N-1)
      complex H(0:N-1)
      Integer kl,kn
      integer hi, hj, i, index(0:n-1), ix

c     initialization of index

!HPF$ INDEPENDENT, LOCAL_ACCESS
      do i = 0, n-1 
         index(i) = i
      end do

c     executing bit reverse for all elements of index

      do L=1,k
        kl = 2**(L-1)       
        kn = N / (kl + kl)    

!HPF$ INDEPENDENT, LOCAL_ACCESS
        do i = 0, n-1 
           ix = index(i)
c          hi = ix / 2
           hi = ishft (ix,-1)
c          hj = hi / kn
           hj = ishft (hi,L-k)
c          hi = mod (hi,kn)
           hi = iand (hi,kn-1)
           index (i) = hi + 2*kn*hj + mod(ix,2)*kn
        end do
       end do
       index = index + 1
       H = C(index-1)
c      call global_get (H, C, index-1)
!HPF$ INDEPENDENT, LOCAL_ACCESS
       do i = 0, n-1 
          C(i) = H(i)
       end do
       end 

C     ********************************************************


C     -------------------------------------------------------
C     ----- Berechnung der Potenzen der Einheitswurzeln -----
C     ----- der Alg. spiegelt eine Rekurrenz erster     -----
C     ----- Ordnung wieder.                             -----
C     -------------------------------------------------------

      subroutine Roots (W,N2,N,Direkt)

      complex W(0:N2-1),omega
cmf$  layout W(:serial)
      real phi
      real Direkt

      phi = Direkt * 8.0 * Atan(1.0) / real(n)
      omega = cmplx(cos(phi),sin(phi))
     
      W(0) = 1
      do i = 1,N2-1
         W(i) = W(i-1) * omega
      end do
      end 



C     ********************************************************


C    **********************************************************


C    --------------------------------------------
C    ----- Transformation nach Cooley-Tukey -----
C    --------------------------------------------

      subroutine Transform(C,W,N,N2,k)

C     Feld c repraesentiert die Fourierkoeffizienten
C     Feld w beinhaltet die Potenzen der Einheitswurzeln

      complex C(0:N-1),W(0:N2-1),H(0:N-1)
cmf$  layout W(:serial)
      complex t

      integer kl,kn,kl1,hi,index(0:N-1),ni
      
      do L=1,k
        kl = 2**(L-1)
        kn = N2/kl
        kl1 = 2 * kl
!HPF$ INDEPENDENT, LOCAL_ACCESS
        do i = 0, n-1
           index(i) = ieor (i,kl)
        end do
        index = index + 1
        H = C(index-1)
!HPF$ INDEPENDENT, LOCAL_ACCESS
        do i = 0, n-1 
c          hi = mod (i,kl1)
           hi = iand (i,kl1-1)
c          ni = mod (i,kl)
           ni = iand (i,kl-1)
           if (hi .ge. kl) then
              C(i) = H(i) - C(i) * w(ni*kn)
            else
              C(i) = C(i) + w(ni*kn) * H(i)
           end if
        end do
      end do
      end 
